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Granular gaseous mixtures under rapid flow conditions are usually modeled as a multicomponent 
system of smooth inelastic hard disks (two dimensions) or spheres (three dimensions) with constant 
coefficients of normal restitution Uij . In the low density regime an adequate framework is provided 
by the set of coupled inelastic Boltzmann equations. Due to the intricacy of the inelastic Boltzmann 
collision operator, in this paper we propose a simpler model of elastic hard disks or spheres subject 
to the action of an effective drag force, which mimics the effect of dissipation present in the original 
granular gas. For each collision term ij, the model has two parameters: a dimensionless factor f3ij 
modifying the collision rate of the elastic hard spheres, and the drag coefficient Both parameters 
are determined by requiring that the model reproduces the collisional transfers of momentum and 
energy of the true inelastic Boltzmann operator, yielding flij — (I + atj)/2 and oc 1 — a?j, where 
the proportionality constant is a function of the partial densities, velocities, and temperatures of 
species i and j. The Navier-Stokes transport coefficients for a binary mixture are obtained from 
the model by application of the Chapman-Enskog method. The three coefficients associated with 
the mass flux are the same as those obtained from the inelastic Boltzmann equation, while the 
remaining four transport coefficients show a general good agreement, especially in the case of the 
thermal conductivity. The discrepancies between both descriptions are seen to be similar to those 
found for monocomponent gases. Finally, the approximate decomposition of the inelastic Boltzmann 
collision operator is exploited to construct a model kinetic equation for granular mixtures as a direct 
extension of a known kinetic model for elastic collisions. 

PACS numbers: 45.70.Mg, 05.20.Dd, 05.60.-k, 51.10.+y 



I. INTRODUCTION 

Natural and industrial granular media are generally present in polydisperse form. In some cases a certain degree of 
polydispersity in masses and/or sizes is unavoidable, while in other cases one is dealing with a real mixture constituted 
by grains belonging to species characterized by distinct mechanical parameters. In conditions of rapid flow, inelastic 
binary collisions are the primary mechanisms and so a kinetic theory description applied to inelastic hard spheres 
has proven to be adequate [1, 2]. In the low density regime, all the relevant information on the state of the mixture 
is contained in the one-particle velocity distribution functions /i(r,v;f,), which obey a set of coupled Boltzmann 
equations, 



N 

(a i +v.V)/ ! =^Jif 1 [v|/,,/ J ], (1.1) 



where N is the number of species and Jfj Cl [v\fi, fj] denotes the inelastic Boltzmann operator that gives the rate 
of change of fi due to collisions with particles of species j. This collision operator depends parametrically on the 
coefficient of normal restitution < 1 (here assumed to be constant) for collisions between particles of species i and 
j- 

Obviously, the problem posed by Eq. (1.1) is much more complicated than in the case of a single granular gas. Not 
only one has to deal with N coupled equations, but in addition the space of parameters is much larger: there are 
N — 1 independent mole fractions, N — 1 mass ratios, N—l size ratios, and N(N — l)/2 coefficients of restitution. An 



important consequence of inelasticity is the breakdown of energy equipartition, even in homogeneous and isotropic 
states. This means that one can associate a different granular temperature to each species [3-5], as confirmed by 
computer simulations [5-9] and real experiments [9, 10] of agitated mixtures. In the case of small spatial gradients, 
the set of Boltzmann equations (1.1) can be solved by means of the Chapman-Enskog method to Navier-Stokes (NS) 
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order. Many attempts to determine the NS transport coefficients are restricted to the quasielastic limit (o^ « 1), 
assuming an expansion around Maxwcllians at the same temperature [11]. A more general derivation takes into 
account the nonequipartition of energy and determines the transport coefficients without any a priori limitation on 
the degree of dissipation [12]. The accuracy of this latter approach has been confirmed by computer simulations in 
the cases of the diffusion [13, 14] and shear viscosity [15] coefficients. 

However, most of the inhomogeneous situations are characterized by a coupling between inelasticity and spatial 
gradients, and so they require a description beyond the NS domain. A typical example of this coupling is represented 
by the simple shear flow [1, 16], where non-Newtonian effects are unavoidable in the steady state. Needless to say, 
the analysis of this type of more general situations based on the inelastic Boltzmann equation (1.1) becomes very 
intricate, especially for multicomponent systems. In order to overcome these difficulties, a possible strategy consists 
of replacing the true Boltzmann collision operator J-" el by a simpler model term that retains its physically relevant 
properties. Following an idea previously proposed for monodisperse granular gases [17], here we explore the possibility 
of describing the multicomponent gas of inelastic hard spheres by a model of elastic hard spheres subject to the action 
of an effective drag force with a different drag coefficient for each species. The parameters of the model are explicitly 
determined by optimizing the agreement with the collisional transfer of momentum and energy obtained from the 
original operator . The resulting model is simpler than the original Boltzmann equation since all the dependence 
on the coefficients of restitution appears explicitly outside the collision term. 

The paper is organized as follows. In Section II the model of elastic hard spheres subject to a drag force is 
formulated, some technical details being relegated to Appendix A. In order to assess the reliability of the model, in 
Section III we present its Chapman-Enskog solution, the expressions of the NS transport coefficients being given in 
Appendix B. In addition, the dependence of the transport coefficients on dissipation is compared with known results 
derived from the Boltzmann equation for inelastic hard spheres [12, 18]. The paper is closed by a discussion of the 
results in Section IV, where a kinetic model equation for granular mixtures is proposed in Appendix C. 

II. PROPOSAL OF THE MODEL 

A. The inelastic Boltzmann equation 

Consider an TV-component mixture composed by smooth inelastic disks (d = 2) or spheres (d = 3) of masses 
and diameters a%, i = 1, . . . , N. The inelasticity of collisions between a sphere of species i and a sphere of species j is 
characterized by a constant coefficient of restitution < a>ij < 1. In the low density regime, the distribution functions 
fi(r, v; t) are determined from the set of nonlinear Boltzmann equations (1.1), where the Boltzmann collision operator 
is 

■4 nel N/i>£] = 4 rl J dv 2 / dse(<7- gl2 )(<7- gl2 ) 

x [ar. 2 / i (r,vi,i)/ J -(r,v 2 ,t)-/ i (r,v 1 ,i)/ J -(r,v 2 ,t)] . (2.1) 

In Eq. (2.1), d is the dimensionality of the system, cry = (oi + aj) /2, <r is a unit vector along the line of centers, is 
the Heaviside step function, and gi 2 = vi — v 2 is the relative velocity. The primes on the velocities denote the initial 
values {v' 1; v 2 } that lead to {vi, v 2 } following a binary (restituting) collision: 

vi = vi - iiji (l + "ij 1 ) (ct • gi2)?, (2.2a) 

v 2 = v 2 + mi (1 + o^ 1 ) {a • gia)?, (2.2b) 

where = m,/ (m, + m,j), so that fiij + fiji = 1. 

The relevant hydrodynamic fields are the number densities rii, the flow velocity u, and the temperature T. They 
are defined in terms of moments of the distributions fi as 

ni= [ dv/i(v), (2.3) 



N N ,. 

pu = J^m^u, = ^mj / dvv/;(v), (2.4) 

i=l i=l 
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N 2 

nT = p = mTi = £ ^ J dv V 2 .h(v), (2.5) 



i=l i=l 



where V = v — u is the peculiar velocity, n = n « i s the total number density, p = X^iLi Pi = Sill TO i n i is the 

total mass density, and p is the pressure. Furthermore, the second equality of Eq. (2.4) and the third equality of Eq. 
(2.5) define the flow velocity and the kinetic temperature Tj for each species, respectively. 
A collision ij conserves the particle number of each species and the total momentum: 



/ 



dvJjf'M/*, /;]=(), (2.6) 



™, J dvvJjfM/,, fj] + mj J dvvJ^Mfj, fi] = 0. (2.7) 
However, unless a^- = 1, the collision ij does not conserve the kinetic energy, so that 

™, J dvv 2 J™ l [v\f u /,] + m j | dv^jjr'M/,-,/,] = -fty, (2.8) 
where > 0. The total "cooling rate" due to inelastic collisions among all species is given by 

so that the rate of change of the granular temperature T due to all the collisions is 



dT 
~dt 



coll 



= -(T. (2.10) 



B. Model of elastic hard spheres with a drag force 

The dependence of the collision operator J]" cl on ctij is rather involved since it appears as the factor a7 2 in the gain 
term and also through the scattering rules (2.2). This represents an additional difficulty with respect to the elastic 
operator Jfj. In order to simplify the cnj -dependence of the inelastic collision operator, we propose the following 
model of elastic particles subject to a drag force [17]: 

JlFMfi, fj] - Pa J$W\fi, fj] + ^ • [(v - Ui) /i(v)] = KaWfi, fr], (2.11) 

where (3ij and Qj arc determined by optimizing the agreement between the model and the true operator. The 
dimcnsionless factor modifies the collision rate of the elastic spheres to mimic that of the inelastic spheres. The 
quantity Qj > is the coefficient of the drag force = — (mj(jj/2) (v — Uj) felt by the elastic spheres of species i. 
This non-conservative force intends to mimic the loss of energy that the true inelastic spheres of species i suffer when 
colliding with spheres of species j. For simplicity, the drag force has been chosen proportional to the velocity 
relative to the mean flow velocity of species i. According to the model (2.11), the Boltzmann equation (1.1) becomes 

1 N d N 

d t fi + v • V/, + — £ »z ■ ( F ijfi) = £/MiM/i. fA- ( 2 - 12 ) 

^ 3 = 1 ^ 3 = 1 

In this way, the problem of a mixture of inelastic hard spheres is mapped, via a renormalization of the collision rate 
and the introduction of a drag force, onto the problem of a mixture of elastic hard spheres. 

The model (2.11) trivially satisfies the mass conservation equation (2.6). In addition, if we assume the symmetry 
relation [3ij = (3ji (to be confirmed later), the momentum conservation equation (2.7) is also verified. Finally, Eq. 
(2.8) yields 



W'iCij'-^i ^jCji^j — <^^^' (2.13) 
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where we have introduced the quantity 



= ri -^(u,-u) 2 . (2.14) 



d 

Note that, according to Eqs. (2.9) and (2.13), the cooling rate ( of the mixture can be expressed in terms of the drag 
coefficients Qj as 



1 N 

C=yE XidjTi, (2.15) 



where /n is the mole fraction of species i. 

Thus far, the parameters fiij and Cij of the model remain unknown, except for the constraint (2.13). In order to 
determine them, we impose that the collisional transfer of momentum and energy of species i due to collisions with 
particles of species j must be the same as those given by the true Boltzmann equation, 

J dvvJlfMfufj] = I dvvKvMfiJj], (2.16) 
J dv^jifM/,,/,] = J dw^^v^Jj]. (2.17) 



This gives 



N ~ /<ivv./;; v. /;.//• (2 ' 18) 



C« = ^ {&J / dv^J^vl/,,/,] - I dv^jjfVl/i,/,]} . (2.19) 
The collision integrals can be simplified by using the property [2] 

J dvMviViflvilfiJj] = 4 rl J dv i / dv 2 /i(v 1 )/ j (v 2 ) 

x J da Q{a • g 12 )(£ ■ g 12 ) [h(v'{) - /i(vi)] , (2.20) 

where /i(v) is an arbitrary test function and 

v'/ = Vl - ^(1 + a«)(? • gi2)5 (2.21) 
is the post-collisional velocity of a particle of species i. Taking /i(v) = v in Eq. (2.20), Eq. (2.18) simply becomes 

N = l ~ ± ^ L , (2-22) 

in agreement with the symmetry relation /3y = (3ji. Next, taking h{y) = v 2 and making use of 

v" 2 -v\ = + aij)(B ■ gi 2 ) [2(ct • G 12 ) + /^(l - ajj)(? • g 12 )] , (2.23) 

where Gi 2 = + /UjjV 2 is the center of mass velocity, one gets from Eq. (2.19) 

<" - nw+W - "«> % " <»*■>«• <2 - 24 > 

Here. 

^«)« = / dv i / dv 2 / J (v 1 )/,(v 2 )^ 2 (2.25) 
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and use has been made of the result [19] 

da-Q(a ■ gi 2 )(5 ; • g i2 ) 3 



■(1 - a^)mn 



,(d-l)/2 



Inserting Eq. (2.24) into Eq. (2.13) we get 



r((d + 3)/2) 



r((d + 3)/2 



mi + mj 



l (.9 12 >, 



(2.26) 



(2.27) 



According to Eq. (2.24), the drag coefficient Qj is a positive definite quantity which only vanishes if the ij collisions 
are elastic (a, 3 = 1). It can be interpreted as the cooling rate of species i due to the inelasticity of collisions with 
particles of species j. To make this more explicit, let us consider the contribution of collisions with particles of species 
j to the rate of change of the partial temperature Tj, i.e., 



dt 



coll,j 



m, 
dn 



"/dv(v u,r./;;' ; v /,../' 



= |8o~ y dv(v - u,) 2 J^ [v|/i, /,] - CyTi, 



(2.28) 



where in the last step use has been made of Eqs. (2.16) and (2.17). The first term on the right-hand side represents 
the contribution to the rate of change not directly associated with inelasticity (except for the presence of the factor 
fy). This term can be either positive or negative, depending on fi and fj. In particular, if = Uj and /j and fj 
are approximated by Maxwellians, the sign of this term is the same as that of the temperature difference Tj — Ti. 
Therefore, the second term on the right-hand side of (2.28) is the genuine contribution associated with inelasticity. It 
must be emphasized that Eq. (2.28), with fyj and Qj given by Eqs. (2.22) and (2.24), respectively, is exact and so it 
is not restricted to the model (2.11). 

While Eq. (2.24) is formally exact, it involves the average {g^ij, which is a functional of fi and fj. An estimate 
of this average can be obtained by assuming Gaussian forms for fi and fj given by 



/i(v) 



v) = n, ( 



\2ttT, 



d/2 



exp 



TOj(v - Ui) 2 



2Ti 



(2.29) 



with an analogous form for fj(v). The distribution (2.29) is the one that, sharing with the exact distribution the 
first d + 2 velocity moments, maximizes the missing information defined as — J dv/j(v) ln/j(v). The corresponding 
average (gfyij 1S obtained in Appendix A by neglecting terms of order fourth and higher in the difference Ui — u 
Inserting Eq. (A7) into Eq. (2.24), one obtains 



3 * 



, _1, ,2 



VtliTj 3 TOj 2 



mjT, 2d T , 



(1 - 0?ij\ 



where 



47r ( d -l)/2 / 2Ti 2T 



1/2 



dT(d/2) 



m, m 7 - 



(2.30) 



(2.31) 



is an effective collision frequency of species i due to collisions with particles of species j. Note that n^j 
while miTiiTidj = mjrijTjQi. Equation (2.30) reduces to the one derived in Ref. 17 when = Uj. 

In summary, the model is defined by the replacement (2.11) with the parameters (3ij and Qj given by Eqs. (2.22) 
and (2.30), respectively. While (3ij depends only on the coefficient of restitution ay, the cooling rate Qj is also a 
function of the masses and sizes of particles of species i and j, as well as of the first few velocity moments (density, 
flow velocity, and partial temperature) of both distributions /, and fj. The expressions of faj and Qj preserve the 
collisional transfer of momentum and energy of the true inelastic Boltzmann equation, although in the latter case we 
have considered the Gaussian forms (2.29) for fi and fj parametrized by their first d + 2 moments in order to get 
explicit results. 

It is worth mentioning that, while the model (2.11) is mathematically simpler than the original Boltzmann operator, 
its functional dependence on fi and fj through the corresponding flow velocities and partial temperatures is highly 
nonlinear. Therefore, the bilinear property 



incl j 



V\fi,fj} 



(2.32) 
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is satisfied, but not the other bilinear property 

JlfMf* +fi2,fji +f j2 ] = JlfMfnJji] +4f 1 [v|/ i i,/ J - 2 ] +4 nel M/ i2 ,/ jl ] + Jlf[v\f i2 ,f j2 ]. (2.33) 

A consequence of the failure of to satisfy the property (2.33) is that, in general, one does not get a closed equation 
for the total distribution function / = J2i h hi the case of mechanically equivalent particles, unless Ui — u and 
Tj = T. However, this drawback is not relevant in most of the situations of physical interest, such as nonequilibrium 
steady states, since in those cases the existence of different velocities and/or temperatures is a consequence of the 
particles being mechanically different. 



C. Homogeneous cooling state 

The simplest application of the model corresponds to the so-called homogeneous cooling state (HCS) [4]. This state 
is characterized by the absence of gradients, so that Ui = u. In that case, the model (2.11) yields 

d t T t = -QT U (2.34) 

where 

& = E J dvV2j $ [vlfi > fj] + c y ■ (2 - 35) 

As said above in connection with Eq. (2.28), Eq. (2.35) is also valid for the inelastic Boltzmann equation, provided 
that Qj is given by Eq. (2.24). 

To evaluate the collision integral in Eq. (2.35), we can take again Eq. (2.29), which (since = u) becomes 



U -> Am = n t (^A ' cxp(-m 4 y 2 /2T J ). (2.36) 



This yields 



Therefore, 



raj 
driiTi 



J ^JllWuM = ^ ^Z/^ - (2 ' 37) 



N - 

_ ^ m t m 3 Tj - Tj I - a tJ ( rrij Tj \ 

Q -h ^ K + ( + lj) [~TT + ~2~ Wi + tJ ■ (2 ' 38) 



i=i 

This expression coincides with the one obtained from the original Boltzmann equation when the approximation (2.36) 
is used. 

The HCS condition is d t Ti/T = 0, which implies d = ( 2 = ■ • ■ = Cn- This gives the N — 1 temperature ratios 
li = Ti/T as functions of the mole fractions, the mass ratios, the size ratios, and the coefficients of restitution of the 
mixture. Comparison with computer simulations [5-7] shows an excellent agreement with the results obtained from 
Eq. (2.38), even for strong dissipation. 



III. NAVIER— STOKES TRANSPORT COEFFICIENTS OF THE MODEL 

To assess the reliability of the model described by Eqs. (2.11), (2.22), and (2.30) in inhomogeneous situations, 
in this Section we will consider the corresponding expressions for the NS transport coefficients of a binary mixture 
(N = 2) and will compare them with the results derived from the original Boltzmann equation [12, 18]. 

The NS coefficients are defined through the constitutive equations [12] 

= _«n w ^ _ l D , VT _ P h = _j (3 1} 

P T p 

Pu = pSke - V ( Vfc"£ + V £ u fe - ^4<?V • u J , (3.2) 
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FIG. 1: Plot of the reduced shear viscosity coefficient ri(a)/ri(l) as a function of the coefficient of restitution a for a three- 
dimensional equimolar mixture in the cases an = a\2 = Q22 = a (case A) and an = a, Q12 = (1 + a)/2, 022 = (3 + a)/4 (case 
B). The panels correspond, from top to bottom, to the systems {m-i/mi, 01/02} = {1,1}, {1,2}, {2, 1}, and {2, 2}, respectively. 
The solid lines are the Boltzmann results, while the dashed lines are the predictions of the model (2.11). 



-T 2 D"\7 Xl - AVT - LVp, 



(3.3) 



where ji = mini (ui — u) is the mass flux of species 1, Pki is the pressure tensor, and q is the heat flux. The transport 
coefficients in the constitutive equations are 



/ D \ 

D' 

D p 
V 

D" 
A 

V L ) 



/ diffusion coefficient \ 
thermal diffusion coefficient 
pressure diffusion coefficient 
shear viscosity 
Dufour coefficient 
thermal conductivity 
\ pressure energy coefficient / 



(3.4) 
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FIG. 2: Plot of the reduced Dufour coefficient D" (a) / D" (f) as a function of the coefficient of restitution a for a three- 
dimensional equimolar mixture in the cases an — at\2 = a.22 = ol (case A) and an = a, Q12 = (1 + a)/2, 022 = (3 + a)/4 (case 
B). The panels correspond, from top to bottom, to the systems \m\jm2, 0"i/o"2} = {1, 2}, {2, 1}, and {2, 2}, respectively. The 
solid lines are the Boltzmann results, while the dashed lines are the predictions of the model (2. If). 



Explicit expressions for the above coefficients are obtained by solving the inelastic Boltzmann equation by means of 
the Chapman-Enskog method. These coefficients are formally given in terms of the solutions of coupled linear integral 
equations involving the linearized Boltzmann collision operators 

= -JirHA (0) ,/i (1) ] - /i (0) ] - ^rU (1 \/ 2 (0) L (3.5) 

M x f¥ = -jg\ff\f ( i-\ (3.6) 

The corresponding expressions for the operators £2 and M2 can be easily obtained from Eqs. (3.5) and (3.6) by just 
making the changes 1 2. In the above equations, is the local version of the distribution function of species i 
in the HCS [4]. 

In the case of the model (2.11), one can follow the same formal steps as in the case of the true Boltzmann equation, 
except that the operators C\ and M.\ are now 

AA (1) = -^(j 1 1[/}°\/ 1 (1 Vj 1 1[/ 1 ( ^/ 1 (0) ])-^^[A ( ^/f] 

.v/W-il-A/W), (3.7) 
2 \ov pi av J 

MifP = -±^JlM \tf\ (3-8) 
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FIG. 3: Plot of the reduced thermal conductivity coefficient A(a)/A(l) as a function of the coefficient of restitution a for a 
three-dimensional equimolar mixture in the cases an = aii = a.22 = a (case A) and an = a, ot\2 = (l + a)/2, a.22 = (3 + a)/4 
(case B). The panels correspond, from top to bottom, to the systems {mi/m2, ^1/02} = {1 , 1}, {1,2}, {2,1}, and {2,2}, 
respectively. The solid lines are the Boltzmann results, while the dashed lines are the predictions of the model (2.11). 

In Eq. (3.7), the drag coefficients (2.30) are those of the HCS, namely 

& = ^"?'( 1+ ^) (1 - 4) ' <M) 

The integral equations defining the transport coefficients are usually solved by expanding in Soninc polynomials. 
For practical purposes, only the leading terms are retained. In addition, we will use the Maxwellian distribution 
(2.36) as a trial function for f^°\ The final expressions for the seven transport coefficients of the mixture are given 
in Appendix B. A symbolic code providing the transport properties under arbitrary values of composition, masses, 
sizes, and coefficients of restitution can be downloaded from the website given in Ref. [20]. 

Let us compare the predictions of the model (2.11) for the NS transport coefficients with the results derived from 
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FIG. 4: Plot of the reduced pressure energy coefficient L(a)n/X(l) as a function of the coefficient of restitution a for a three- 
dimensional equimolar mixture in the cases an = «i2 = Q22 = a (case A) and an = a, ai2 = (l + a)/2, a22 = (3 + a)/4 (case 
B). The panels correspond, from top to bottom, to the systems {mi/m2,ffi / '02} = {1,1}, {1,2}, {2, 1}, and {2, 2}, respectively. 
The solid lines are the Boltzmann results, while the dashed lines are the predictions of the model (2.11). 



the original inelastic Boltzmann equation [12, 18, 21]. The three coefficients D, D', and D p associated with the mass 
flux are identical in both descriptions. This is a consequence of the requirement (2.16). For this reason, these three 
coefficients will not be shown here. Their behaviors for several representative cases have been analyzed elsewhere 
[18]. The other four coefficients differ in both descriptions, except trivially in the elastic case (ay = 1). Since 
the parameter space of the problem is six-dimensional, namely {xi, mi/m2, cri/172, an, «i2, CV22}, it is convenient to 
choose some specific cases. First, we consider hard spheres (d = 3) with a common coefficient of restitution, i.e., 
an = ai2 = 022 = a, (case A) and also with the choice of the coefficients of restitution an = a, ai2 = (1 + a)/2, 
a22 = (3 + a)/4 (case B). Note that in case B one has an < ai2 < 022- Consequently, for given values of x\, mi/m2, 
C1/CT2, and an = a, the system A is more inelastic than the sytem B. Next, we restrict ourselves to equimolar 
mixtures {x\ = ^). This reduces the parameter space to three quantities, namely {mi/m2, <J\l<Ji, a}. To focus on 
the influence of inelasticity on the transport coefficients, we fix the values of the mass and size ratios and plot the 
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transport coefficients as functions of a. In addition, the coefficients are reduced with respect to their values in the 
clastic limit, except in the case of the pressure energy coefficient L, which vanishes for elastic collisions if mj = m 2 . 
In this case, the plotted quantity is L(a)n/X(l), where A(l) denotes the elastic value of the thermal conductivity 
coefficient. As representative cases we have chosen {™i/m 2 , ai/a 2 } = {1, 1}, {1, 2}, {2, 1}, and {2, 2}. Note that the 
system {mi/m 2 , 01/0-2} = {1,1} corresponds to mechanically equivalent particles in case A. 

The shear viscosity is shown in Fig. 1. We observe that, at least for the cases analyzed here, the model underestimates 
the Boltzmann values. Although the relative discrepancies increase with dissipation, they are practically insensitive 
to the mass and size ratios. As a consequence, the model captures well the influence of mi/m 2 and 01/02 on r\. In 
particular, it is interesting to note that the ratio rj(a)/r](l) for the system {mi/m 2l 01/02} = {1,1} is practically 
indistinguishable from that of the system {mi/m 2 , 01/02} = {2,2}, this effect being reproduced by the model. As 
expected, the influence of dissipation is less significant in case B than in case A. 

Now we consider the three coefficients associated with the heat flux. We start with the Dufour coefficient D" , 
which is plotted in Fig. 2. The system {mi/ra2, 01/02} = {1, 1} is not displayed since D" = at any value of a for 
mechanically equivalent particles (case A). In case B, however, D" ^ for the system {toi/to2, 01/02} = {1, 1} since 
an ^ ol\2 a 22 . We have checked in that case that the performance of the model is quite good. It is apparent from 
Fig. 2 that, for a given value of a, the ratio D"(a)/D"(l) has a significant dependence on mi/m 2 and/or 01/02. 
This effect is well accounted for by the model. It is noteworthy the dramatic influence of inelasticity on the value of 
the Dufour coefficient when mi — m 2 and 01 ^ a 2 . As shown in the top panel of Fig. 2, this feature is accurately 
predicted by the model. 

The thermal conductivity and the pressure energy coefficients are displayed in Figs. 3 and 4, respectively. In contrast 
to D", these coefficients are meaningful in the case of mechanically equivalent particles. The model performs quite 
good a job, even for strong dissipation, except perhaps for the most disparate mixture {mi/m 2 , 01/02} = {2,2} in 
case A. 

In summary, it is fair to say that the model of elastic spheres subject to a drag force mimics, at least at a semi- 
quantitative level, the influence of inelasticity on the transport coefficients of a binary mixture of inelastic hard 
spheres. More specifically, the mass flux coefficients (D, D 1 , and D p ) are the same in both systems, while the shear 
viscosity 77 is underestimated. In the least favorable case (case A), the discrepancies of the values of r\ are about 7% 
at a — 0.8 and 13% at a — 0.6 for the systems studied here. With respect to the heat flux coefficients (D" , A, and 
L), the general agreement is, paradoxically, better than in the case of the shear viscosity. For instance, the deviations 
of the thermal conductivity are less than 3% at a — 0.8 and 8% at a — 0.6 for the systems studied here in case A. 
It is also interesting to remark that the reliability of the transport coefficients predicted by the model is practically 
independent of the disparity in mass and size ratio. 

IV. DISCUSSION 

The study of multi-component granular systems is of paramount importance from a practical point of view, but 
also at a fundamental level. In the low density regime, granular gases are well described by the inelastic Boltzmann 
equations with constant coefficients of restitution < 1. However, the intricacy of the inelastic Boltzmann collision 
operator J™ ' makes it difficult to extract explicit results, especially at moderate or strong dissipation. This motivates 
the search for models that, while retaining the basic properties of Jlf , provide a simpler framework for granular 
mixtures. 

In this paper, we have proposed a model of elastic hard spheres under the action of an external drag force. The 
role of the force is to mimic the collisional energy loss in the true gas of inelastic spheres. In addition, the collision 
rate for elastic collisions is modified by a factor /3jj with respect to the one for inelastic collisions. Accordingly, in 
this model the collision operator J-" el is replaced by the operator Kij defined in Eq. (2.11). While the dependence 
of J-" cl on OLij appears in the gain term and through the collision rules [see Eqs. (2.1) and (2.2)], such a dependence 
appears explicitly in through the parameters /3,j and Qj . These parameters are determined by imposing that the 
model operator reproduces the correct collisional transfer of momentum and energy, Eqs. (2.16) and (2.17). The 
former condition gives the simple expression (2.22) for the parameter /3ij modifying the collision rate. The energy 
condition (2.17) shows that the drag coefficient Qj is a functional of the distribution functions fi and fj that is simply 
proportional to 1 — cqj, Eq. (2.24). Thus, Qj can be interpreted as the cooling rate of species i due to collisions with 
particles of species j, as indicated by Eq. (2.28). Since Qj requires the average {g\ 2 )ij defined by Eq. (2.25), the 
Gaussian forms (2.29) parametrized by the partial velocities and temperatures are used to estimate Qj, yielding Eq. 
(2.30). 

By construction, the model leads to the same results for the temperature ratios in homogeneous states as those 
obtained from the inelastic Boltzmann equation in the multi-temperature Maxwellian approximation [4, 5], which 
is known to compare quite well with computer simulations [5-7]. To measure the performance of the model in 
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inhomogeneous situations, the seven NS transport coefficients for a binary mixture predicted by the model have been 
computed and compared with previous results from the inelastic Boltzmann equation [12, 18, 21]. The three transport 
coefficients associated with the mass flux are identical in both descriptions in the first Sonine approximation, while 
the remaining four coefficients differ, as illustrated in Figs. 1-4. We found that the model captures reasonably well the 
dependence of the shear viscosity r\ and the three heat flux coefficients D", A, and L on dissipation. This agreement 
is especially remarkable in the case of the thermal conductivity A. The degree of agreement between the model and 
the inelastic Boltzmann transport coefficients has been found to be practically independent of the disparity of mass 
and/or size, being similar to the one found for the one-component case [17]. Beyond the NS domain, where the spatial 
gradients are not sufficiently small, it is uncertain whether or not the model predictions are close to the inelastic 
Boltzmann ones. However, at least in the uniform shear flow problem, the simulation results in the one-component 
case are practically indistinguishable in both systems [22]. We expect that such a good agreement is also kept for 
multicomponent systems, given that the discrepancies between the model and the inelastic Boltzmann transport 
coefficients are not significantly affected by the mass and size ratios. 

One of the main features of the model (2.11) is that it allows, in an approximate way, to "disentangle" in Jf- el the 
purely dissipative effects (represented by the drag- force term) from those (represented by PijJfj) which are essentially 
present in the elastic case. Taking advantage of this decomposition, one can extend to the inelastic case any model 
kinetic equation proposed for ordinary multicomponent gases. This is quite an important issue since, even for elastic 
collisions, the Boltzmann equation is too complex to study far from equilibrium situations [23], for which the NS 
description fails. For one-component ordinary gases, the prototype model kinetic equation is the Bhatnagar-Gross- 
Krook (BGK) model [24, 25]. Several extensions of this model to the case of one-component inelastic hard spheres 
have been proposed [26-28]. Furthermore, some kinetic models for ordinary gas mixtures inspired on the BGK model 
can be found in the literature [23, 29-31]. The common structure of the latter models is 

^[v|/i,/i]^-^[/i(v)-&(v)], (4.1) 

where Vij is a velocity-independent effective collision frequency of a particle of species i with particles of species 
j and /ij(v) is a reference distribution function whose velocity dependence is explicit and that involves a number 
of parameters to be determined by imposing that Eq. (4.1) retains the main physical properties of the Boltzmann 
operator Jfj. Using the mapping (2.11), any kinetic model of the form (4.1) can be easily extended to the case of 
inelastic collisions as 

4 n< Vl/i, fj] - -PVii [/i(v) - fij(v)] + ^ • [(v - Ui) Mv)] , (4.2) 

where /3y and Qj are given by Eqs. (2.22) and (2.30), respectively. A specific choice for based on the Gross-Krook 
kinetic model [29] is worked out in Appendix C. Several applications of the kinetic model (4.2) to non-Newtonian 
flows are under way and will be published elsewhere. 

Acknowledgments 

This research has been supported by the Ministerio de Education y Ciencia (Spain) through Programa Juan de la 
Cierva (F.V.R.) and Grant No. FIS2004-01399, partially financed by FEDER funds. 

APPENDIX A: EVALUATION OF ( S ? 2 ) i ,-, (siagu)^, AND (jugu, • Gi 2 >y 

The aim of this Appendix is to evaluate the averages {gi 2 )ij, (<7i2gi2)ij, and (g , i2gi2 ■ G\2)ij by assuming the 
Gaussian forms (2.29). 

Let us start with the evaluation of (gfyij- Inserting Eq. (2.29) into Eq. (2.25), one gets 

<*>«-(af ferwM-^-^)^ + °'- ,,/ ' <ai) 

where wi = vi — u,, w 2 = v 2 — Uj, and W12 = Wi — w 2 . Let us also introduce the variable 
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so that dwidw 2 = dwi 2 dWi 2 and 



miw\ mjW% 
T T 



Ttii m 



T 



3 1 w? 



12 



Ti . To- 



rn,' m,- 



wj 2 . 



(A3) 



Thus, after making the changes (v x ,v 2 ) — ► (w 12 ,W 12 ) and integrating over W 12 , Eq. (Al) becomes 



(9l2) ij 



-d/2 



2T, 2T, 



rrii m 



2T, 2T 7 - \ f , 

1 / dw 12 exp 

'", rrij J J 

'2 

K~ d/2 J dce-^lc + u^-l 3 , 



2Tj 2Ti 



rrii 



-12 



W12 + Uj - U, 



(A4) 



where in the second step we have made the change of variable c = (2Tj/m, + 2Tj/m,j) ' W12 and have defined 
(2Ti/nii + 2Tj/m,j)~ 1 / 2 (ui — Uj). The Gaussian integral in Eq. (A4) cannot be evaluated analytically for 



uj 



arbitrary u*j. Here we will evaluate it by neglecting terms of fourth and higher order in u*^. In that case, one has 

, 2 



: + 14 1 3 = c 3 + 3cc • u?. + |cut- 2 + fc- 1 (c • + • • • , 
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where we have taken into account that the third order terms do not contribute to the integral. As a consequence 

3d + 1 



(A5) 



7r- d / 2 / dce- c |c + u 



* |3 



v- d ' 2 J dce- 

r((d+l)/2)d + l 



3 o u -r i 2 



r(d/2) 



11^*2 
1 + - Uij 



Insertion of (A6) into (A4) yields 



Following similar steps, it is straightforward to get 



1/2 



Ti , Tj 3 2 
rrii rrij M 



(5l2gl2>ij = 



27; + 27} 
m, to. 



dce- c Ic + u^Kc + u^-). 



Next, one has 

|c + u?.|(c + u*j) = c (u*j + c- 2 u* 3 • cc + • • • ) , 
where the ellipsis denote terms which are odd in c or are of at least of third order in u* 3 . Thus, 

d+ 1 



(A6) 



(A7) 



(A8) 



(A9) 



dee c |c + u*|(c + u*) -> u 
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y d 

d + ir((d + i)/2) 
d r(d/2) 



d / 2 / dce- c c 



u 



From Eqs. (A8) and (A9) it follows that 



(5l2gl2)«j = 



d+1 r((d+ l)/2) / 2T, 2Tj 



1/2 



d 



r(d/2) 



(Uj-Uj)- 



(A10) 



(All) 



Finally, we consider the evaluation of (#i 2 gi 2 • Gi 2 )y . Taking into account Eq. (A2) and the definition of Gi 2 [see 
below Eq. (2.23)], it is easy to get the relationship 



rrii rrii 



G 12 =W 12 + + 



rrii 
Ti 



rriirrij 
rrii + rn. 



(Tr^-T^) El2 



(A12) 
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As a consequence, 



(#12gl2 • Gi2>ij = I ^T + ^r 



(ffl2gl2) 4J ' ( + — Uf 



mj 
T- 



mirrij 
rrii + nij 



(A13) 



where, by symmetry, (<?i2gi2 • Wi 2 ),j = when evaluated with the distributions (2.29). Making use of Eqs. (A7) and 
(All) one gets 



(5l2gl2 ' Gl 2 )i 



d+ir((d+l)/2) / 2Tj 2Tj\ 1/2 ( mj mA 
d T{d/2) m/mj { % + % ) 



m, 



Uj + ^U f • (Uj - Uj) 



d — V^-fa - T r ) 

m, + rrin V J / 



-Zjf 3 ^ .2 

1" + ( u * - u j) 

mi TOj 2d 



(A14) 



APPENDIX B: EXPLICIT EXPRESSIONS FOR THE TRANSPORT COEFFICIENTS 

In this Appendix, we take advantage of the results derived from the original Boltzmann equation [12, 21] to 
determine the expressions corresponding to the model (2.11), by using the forms (3.7) and (3.8) for the operators Li 
and Mi. 



1. Mass flux 



In the case of the transport coefficients associated with the mass flux, the results are [12, 21] 

-l 



m,\m 2 n \ 2 



P ^^) p , T + P {^)j D ^ D '\ 



(Bl) 



7i 



(B2) 



D ' = -Yv D v 



(B3) 



Here, ( = d — ( 2 is the cooling rate of the mixture in the HCS [see Eq. (2.38)] and the collision frequency v is given 
by 



The result is 



v = ~T~~^r — / dvV ' ( A/i,mV - ^M 1 f 2 , M V ) ■ 
dn x T~ix J \ 2:272 



v = Vq ^ta\ (1 + Ql2 ) — — F2M21 + Z1/Z12) ■ 



dT(l) 



6162 



(B4) 



(B5) 



where v = na^ 1 \j2T(m\ + m2)/mim 2 , 6\ = (/X2i7i) _1 , #2 = (Mi272) _1 , and we recall that 7, = Tj/T. The 
expressions of the transport coefficients D, D p , and D' are exactly the same as those obtained from the true Boltzmann 
equation. 
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2. Pressure tensor 



The shear viscosity coefficient 77 can be written as 

where the expression of the (dimensioniess) partial contribution 77* is 

7 2 (2t 22 - C*) - 271 T12 



Vi = 2- 



J 7l72[C* - 2 C*(ni + T 22 ) + 4(tiiT 22 - Tl2T2l)] ' 

Here, (* = (/v and we have introduced the (reduced) collision frequencies 

1 



(d-l)(d + 2) mT 2 7 



^-2— / dvR 1M d {h,MRi,ke) 
2 7i^o J 
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where 
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A similar expression can be obtained for fy 2 by just making the changes 1^2. Taking advantage of the results 
derived for the d-dimensional Boltzmann equation [21], one gets for the model the expressions 



Til = 



4^-1)72 



Xi 



d-l 



( — ) (20i)- 1 /Ml + «n) 
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+x 2M2 i(l + ai2 )6l /2 6- 1/2 [(d + 3)( M i 2 e 2 - M21 W 2 (#i + 2)" 1/2 



(Bll) 



4 7r (<*-i)/2 
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X 2 0/ 6 2 ' (1 + «1 2 ) 



d(d + 2)r(f) M12 
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+d^ 21 9- 2 (9 1 + 9 2 f 2 2d{ i^ 1 \: A e 2 \e l + 9 2 r 1/2 



2(d-l) 



where = Cij/^o an d Cj i s given by Eq. (3.9). 



(B12) 



3. Heat flux 

The transport coefficients appearing in the heat flux (D" , L, and A) can be written as [12, 21] 

d + 2 n 



D" 



L = — 



2 (mi + m 2 )vo 



d + 2 T 
2 (mi+m 2 )^o 



M12 M21 VM12 A<21. 
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(B13) 
(B14) 
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d + 2 nT 

A — 



2 (mi + m 2 )^o 
where 



Ml2 M21 \Ml2 M21, 



(B15) 



D = pT D*, D p = —D* D' = —D'*, (B16) 

the coefficients D, D p , and D' being given by Eqs. (B1)-(B3), respectively. The expressions of the (dimensionless) 
coefficients d*, £*, and A* are 

d\ = i{2 [2 Xl2 r 2 - Yi(2 X22 -3C)][xi2X2i-XiiX22 + 2(xii+X22)C-4C 2 ] 

+2 (|£) ^ (Y 3 + Y 5 ) [2 Xl2X2 i + 2xh - C(7X22 - 6C*)] 

-2xi2 (n + r 6 ) (2 X ii + 2x22 -7C*)}, (B17) 

*I = ^ {-2^3 [2(xi2X2i - XHX22)X22 + C*(7xnX22 - 5xi2X2i + 2 X ^ 2 - 6xiiC* - 7X22C* + 6C 2 )] 
+2r 4 xi2 [2x12X21 - 2xiiX22 + 2C(xn + X22) - C 2 ] 

+2Y 5 C [2x12X21 + X2 2 (2X22 - 7C) + 6C* 2 ] - 2xi 2 C*^6 [2(xn + X22) - 7(1} , (B18) 

K = £ {-YsC [2X12X21 + X2 2 (2X22 - 7C) + 6C* 2 ] + Xi2C*n [2(xn + X22) - 7(*} 

-Y B [4 X i2X2i(X22 - C) + 2X2 2 (5C ~ 2 X n) + 2xii(7x2 2 C* - 6C* 2 ) + 5C* 2 (6C* - 7 X 2 2 )] 

+Xi2>6 [4X12X21 + 2xn(5C* - 2x22) + C*(10X22 - 23C*)]} . (B19) 

In the above equations, we have introduced the quantities 

A = [4(xi2X2i - X11X22) + 6C*(xn + X22) - 9C* 2 ] [X12X21 - X11X22 + 2C(xn + X22) - 4(* 2 ] , (B20) 
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The explicit expressions for the coefficients Xij an d w !3 arc 

(d-l)/2 4 / \ d-l 
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(d-l)/2 



+ - 



1 
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02(61 + 2 ) 



3/2 r 
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F + (d + 2 )^±^B 

02 



(B29) 



where 



^ = T (f) /2 d(d + 2) Xl ^ 2l(1 + Ql2)(gl + h)- 1/2 0lV 2 (g^- +Cn +cr 2 , 



A = (d + 2)(20i 2 +0 2 )+4/i 21 (0 1 +0 2 )( ( i-l)^i 2 0r 1 
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(B30) 



(B31) 
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-3(d + S)^^ 1 - /^(d + 3)02 ^fc + 02) 2 + (d + 2)(0! + 2 ), 



(B32) 



Mli0i" 2 (0i + 2 ) 2 (d + 3) [(d + 2)0! + (d + 5)0 2 ] 

+4(rf - 1)^21 (01 + 02 ) {0i20r 2 [(rf + 2)0! + (d + 5)0 2 ] + 20207 1 } 

+3(d + 3)0 2 2 0r 2 [(d + 2)0i + (d + 5)0 2 ] +20 12 1 - 1 [(d + 2) 2 1 + (24 + lid + rf 2 )0 2 ] 

+ (d + 2)0 2 0r 1 [(d + 8)0! + (d + 3)0 2 ] - (d+ 2)(0i + 2 )0r 2 2 [(d + 2)0i + (d + 3)0 2 ], 



(B33) 
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(B34) 



Here, <j>i 2 — ^i 2 2 — A*2i0i- The expressions for X22, X21, an d ^21 can be easily obtained from Eqs. (B28)-(B30) by 
exchanging 1 2. 



APPENDIX C: A MODEL KINETIC EQUATION FOR GRANULAR MIXTURES 

This Appendix addresses the construction of the kinetic model (4.2) for a specific choice of the reference distribution 
function fy. 

The most natural way of extending the BGK model to (elastic) mixtures is to assume a Gaussian form for the 
reference distribution fy, i.e., 

where Ujj and T t j are parameters to be determined. This is the form of the model proposed by Gross and Krook 
[29], which has been widely used in the literature. The usual criteria to determine the unknowns Vij, Ujj, and 



m , 2 



(Cl) 
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are to impose that the kinetic model reproduces the collisional equations of momentum and energy of the original 
Boltzmann operator Jfj, namely, 



j dwJfj[v\fi,fj] = -VijJdwtfiW-fijiv)] 
= -VijUi {Ui - Uij) , 

— (fi - Ty) + 



(C2) 



2 2 

^ - u^ 



Conservation of total momentum in collisions ij implies that 

v lj n l m i (uj - Ujj) + /',,/' ,///, (uj - Uji) = 0. 



(C3) 



(C4) 



For physical reasons it is assumed that rnuij — rijVji, as happens with the effective collision frequencies £y defined 
by Eq. (2.31). Moreover, the reference velocity u i:; - is assumed to be symmetrical, i.e., uy = Ujj. Under the above 
conditions, Eq. (C4) yields 



Uij — jlij W.i -\- [Aji Uj 



(C5) 



To identify the remaining parameters Vij and through Eqs. (C2) and (C3) we need to compute the corresponding 
collision integrals associated with Jfj. In the original model proposed by Gross and Krook [29], the left-hand sides of 
Eqs. (C2) and (C3) were evaluated by considering mixtures of Maxwell molecules, in which case the collision rate is 
independent of the relative velocity. Here, however, we want to keep in Jfj the velocity dependence of the collision 
rate characteristic of hard spheres. 

Applying the property (2.20) and making use of Eqs. (2.21) and (2.23) with atj = 1, it is easy to obtain 
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27r (d-i)/2 
dv 1 v 1 Jf[v 1 \f i ,f j ] = - - 
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lijiCTij ninj(g 12 gi2}ij, 
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J dVivljfj^lfiJj] 
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r((d + 3)/2) 



Hji<Jij ninj(gi 2 gi2 ■ G 12 )i 



(C7) 



where the integrations over the solid angle have been performed. The averages appearing in Eqs. (C6) and (C7) can 
be evaluated, as in the case of (gf 2 )iji by assuming Gaussian forms for fi and fj and neglecting terms of order third 
and higher in — Uj. This is done in Appendix A with the results (All) and (A14), so that Eqs. (C6) and (C7) 
become 



J dvi viJfj[vi\fi, fj] = -ZijriiUjiiui - Uj), 



(C8) 



— + ^ 

rrii rrii 



J dv 1 v^Jfj[v 1 \f i ,f j ] = -liijuaiji 



rrii + rrij 
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rrii rrij 2d 
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respectively, where £y is given by Eq. (2.31). Substitution of Eqs. (C8) and (C9) into Eqs. (C2) and (C3) yields 



Vij — £ij , 



(CIO) 
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Note that — Tj + Tji — Tj = [miirij / (mj + mj)](iij — Uj) /d. It is also interesting to remark that the term 
proportional to (iii — Uj) 2 (Tj — Tj) in Eq. (Cll) is absent in the kinetic model based on the Boltzmann equation for 
Maxwell molecules [23]. 

Equations (C5), (CIO), and (Cll), along with Eq. (CI), close the construction of the kinetic model (4.1) for mixtures 
of elastic hard spheres. Next, the corresponding kinetic model for granular mixtures is defined by Eq. (4.2). This 
extended kinetic model has the same collisional transfer of momentum and energy as the true inelastic term JJ™ , as a 
consequence of Eqs. (2.16), (2.17), (C2), and (C3), at least in the Gaussian approximation (2.29). As a consequence, 
the temperature ratios 7j = Tj/T of the HCS are the same as those obtained from the inelastic Boltzmann equation. 

The NS transport coefficients predicted by the kinetic model can be evaluated from the Chapman-Enskog method. 
Since \\i — u = ji/rmm, one has 



u 



Tii (nii + rrij ) rij {mi + mj ) 



(C12) 



Furthermore, to NS order, T i =T i = 7,T and 



7,+ 



2wim 7 - 



(m; + mj)'- 



{lj - 7») 
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As a consequence, the reference distribution /y becomes 



/«(v)=/ff ) (V) 



(C14) 



where 



/ff } (V)=n s 



27rTi, 



d/2 



exp 



2T, 



-V 



Therefore, in the binary case, the linear operators d and A^i take the forms 
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'-12 



J_l_ 

ni 



A V f (1) -ji — f (0) 



(C16) 



M/'^-^^'^vi (C17) 

^ J 12 ^2 

The transport coefficients are given again by Eqs. (B1)-(B3), (B6), (B7), and (B13)-(B15), but now the associated 
collision frequencies (B4), (B8), (B9), and (B24)-(B26) are computed by using the linear operators (C16) and (C17). 
In particular, the transport coefficients associated with the mass flux are the same as those obtained from the inelastic 
Boltzmann equation [12, 21]. 
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